ohibc logo
OHI British Columbia | OHI Science | Citation policy

knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
                      echo = TRUE, message = FALSE, warning = FALSE)

# library(ggmap)
library(rgdal)
library(sf)
library(broom)

source('~/github/ohibc/src/R/common.R')  ### an OHIBC specific version of common.R

dir_git     <- '~/github/ohibc'
dir_spatial <- file.path(dir_git, 'prep/_spatial')  ### github: general buffer region shapefiles
dir_anx     <- file.path(dir_M, 'git-annex/bcprep')

### goal specific folders and info
goal     <- 'ao'
scenario <- 'v2017'
dir_goal <- file.path(dir_git, 'prep', goal, scenario)
dir_goal_anx <- file.path(dir_anx, goal, scenario)

### provenance tracking
library(provRmd); prov_setup()

### set up the default BC projection to be BC Albers
p4s_bcalb <- c('bcalb' = '+init=epsg:3005')

1 Summary

This data is incorporated into the OHI British Columbia First Nations Resource Access and Opportunities (AO) goal.


2 Data Source

Fishing licenses

Fishery Management Area boundaries

  • Reference: From Karen Hunter, DFO
  • Downloaded: May 24, 2016
  • Description: Pacific Fisheries Sub Management Areas
  • Native data resolution: PFMSA
  • Format: ESRI shapefile

Census subdivision administrative boundaries

2016 Census Total Population Results Census Subdivisions


3 Methods

3.1 Read in data

Read in DFO contamination closures from .xlsx files.

lic_clean_file <- file.path(dir_goal_anx, 'int/licenses_clean.csv')

license_file <- file.path(dir_anx, '_raw_data/fishing_license/Pacific Region Commercial Fishing Licenses En.csv')
metadata_file <- file.path(dir_anx, '_raw_data/fishing_license/Commercial License Meta Data EN (1).csv')

if(!file.exists(lic_clean_file)) {
  
  licenses_raw <- read_csv(license_file) %>%
    setNames(tolower(names(.)) %>% str_replace_all('[\\s()]+', '_'))
  
  
  for(tmp in names(licenses_raw)) {
    ### tmpname <- names(x)[1]
    licenses_raw[[tmp]] <- ifelse(is.na(licenses_raw[[tmp]]) & class(licenses_raw[[tmp]]) == 'character', 
                                  '', 
                                  licenses_raw[[tmp]])
  }
  
  metadata <- read_csv(metadata_file, skip = 54) %>%
    setNames(tolower(names(.)) %>% str_replace_all(' ', '_'))
  
  licenses_clean <- licenses_raw %>%
    left_join(metadata, by = c('licence_prefix' = 'prefix')) %>%
    mutate(lic_holder = ifelse(licence_holder_type == 'VESSEL',
                               paste('vessel', licence_holder_vessel),
                               paste(licence_holder_first_name, licence_holder_midddle_name, licence_holder_last_name)),
           lic_oper   = ifelse(licence_holder_type == 'VESSEL',
                               paste('vessel', licence_operator_vessel),
                               paste(licence_operator_first_name, licence_operator_midddle_name, licence_operator_last_name)),
           vessel_own = ifelse(licence_holder_type == 'VESSEL',
                               paste(vessel_owner_first_name, vessel_owner_midddle_name, vessel_owner_last_name),
                               NA)) %>%
    mutate(year = licence_suffix,
           fn = str_detect(tolower(prefix_description), 'aborig|northern native|northernnative'),
           fn = ifelse(is.na(fn), str_detect(tolower(licence_holder_type), 'aborig'), fn)) %>%
    select(year, 
           lic_prefix = licence_prefix, 
           prefix_desc = prefix_description,
           lic_tab = licence_tab,
           fn, area,
           fishery_group, foreign_vessel,
           lic_hold_type = licence_holder_type,
           lic_hold_vessel = licence_holder_vessel,
           lic_oper_vessel = licence_operator_vessel,
           lic_holder, lic_oper, vessel_own,
           vessel_length_m = vessel_length_mr_)
  
  write_csv(licenses_clean, lic_clean_file)
} else {
  message('File exists: ', lic_clean_file)
  git_prov(c(license_file, metadata_file), 'input')
  git_prov(lic_clean_file, 'output')
}

3.1.1 License areas

This link provides a look-up of different commercial management areas to PFMAs and PFMSAs. This might be useful for spatializing the licenses.

http://www.pac.dfo-mpo.gc.ca/fm-gp/licence-permis/areas-secteurs-eng.html

The table in this link has been saved as a .csv in the ao/v2017/raw folder.

3.2 First Nations license holders by fishery group

This bar chart shows the proportion of licenses allocated to First Nations for each fishery group. The indicated proportion is the maximum for that fishery group across all years.

lic_clean <- read_csv(lic_clean_file)

license_fn_yr <- lic_clean %>%
  select(year,
         prefix_desc, 
         fn,
         lic_prefix, 
         fishery_group) %>%
  mutate(fishery_group = tolower(fishery_group)) %>%
  group_by(year, fishery_group) %>%
  summarize(n_total = n(),
         n_fn = sum(fn),
         pct_fn = n_fn / n_total) %>%
  ungroup()

license_fn_yr <- license_fn_yr %>%
  group_by(fishery_group) %>%
  mutate(max_pct_fn = max(pct_fn, na.rm = TRUE)) %>%
  ungroup()

bar_by_gp <- license_fn_yr %>%
  select(max_pct_fn, fishery_group) %>%
  filter(max_pct_fn > 0) %>%
  distinct() %>%
  arrange(desc(max_pct_fn)) %>%
  mutate(fishery_group = factor(fishery_group, levels = unique(fishery_group))) %>%
  ggplot(aes(x = fishery_group, y = max_pct_fn)) +
    geom_bar(stat = 'identity') +
    theme(axis.text.x = element_text(angle = 75, hjust = 1))

print(bar_by_gp)

3.3 First Nations license ownership

This examines the total fishing license allocation to First Nations across all license types, fishery groups, and regions, by year.

lic_clean <- read_csv(lic_clean_file)

license_fn_yr1 <- lic_clean %>%
  select(year,
         prefix_desc, 
         fn,
          lic_prefix,
          fishery_group) %>%
  mutate(fishery_group = tolower(fishery_group)) %>%
  group_by(year) %>%
  summarize(n_total = n(),
         n_fn = sum(fn),
         pct_fn = n_fn / n_total) %>%
  ungroup()

lic_total_plot <- ggplot(license_fn_yr1, aes(x = year, y = n_total)) +
  geom_area(fill = 'grey60', color = 'grey40', size = .25) +
  geom_area(aes(y = n_fn), fill = 'steelblue3', color = 'grey40', size = .25) +
  labs(y = 'Number of licenses')

lic_pct_plot <- ggplot(license_fn_yr1, aes(x = year, y = pct_fn)) +
  geom_line(color = 'steelblue3', size = .5) +
  labs(y = 'Percent FN licenses')

print(lic_total_plot)

print(lic_pct_plot)

3.4 Spatialize licenses

Much of this code was developed for the AO Closures prep script; since there are likely idiosyncracies within the data here compared to the closures data, the code has been copied and customized rather than functionalized.

3.4.1 Clean areas and subareas

Clean up area and subarea calls from the cleaned data. No need to keep the license info at this point; can rejoin later. Many of the descriptions are complex text such as “Area 12, Area 13, except SA 13-1 to 13-12, 13-15 to 13-17, Area 14 (except SA 14-5, 14-8, 14-15, Area 15 (except SA 15-5), Areas 16 to 29 (except SA 17-20, 23-6 and 29-5)” which must be parsed down to areas.

  • Remove area-subarea notation and just keep the area level
  • In many cases, areas are listed as, e.g., Areas 3 to 10.
    • Expand these to e.g. Areas 3 4 5 6 7 8 9 10.
    • These individual area numbers are then extracted as individual observations.
  • bind with rows where no “to” clause existed
  • join with lookup table of PFMA area to OHIBC rgn
lic_area_all <- read_csv(lic_clean_file) %>%
  select(area, year) %>%
  group_by(area) %>%
  mutate(year_1 = min(year, na.rm = TRUE),
         year_n = max(year, na.rm = TRUE)) %>%
  ungroup() %>%
  select(-year) %>%
  distinct() %>%
  mutate(area = ifelse(is.na(area), 'no area designated', tolower(area)))
  
lic_area_lookup <- read_csv(file.path(dir_goal, 'raw/license_area_lookup.csv')) %>%
  mutate(lic_area = tolower(lic_area),
         desc     = tolower(desc))

# x <- lic_area_all %>%
#   filter(area %in% lic_area_lookup$lic_area) %>%
#   arrange(area)
# y <- lic_area_all %>%
#   filter(!area %in% lic_area_lookup$lic_area) %>%
#   arrange(area)
# z <- lic_area_lookup %>%
#   filter(!lic_area %in% lic_area_all$area) %>%
#   arrange(lic_area)
### All the remaining unmatched license areas (y) seem to have ended long ago,
### most before 2000, and only a couple as late as 2002 or 2003.
### Presumably these license areas are no longer valid, so don't show up
### in the license area designations.

### Because we are counting within regions rather than area-weighting, it is
### less important to worry about subareas.  Drop them for simplicity.  This
### also means not having to worry about "except" or "excludes"
lic_lookup_clean <- lic_area_lookup %>%
  select(lic_area, desc) %>%
  mutate(area_edit = desc,
         area_edit = str_replace_all(area_edit, '(?<=[0-9])[-—.](?=[0-9])', '~'), ### this line writes subareas as ##~##
         area_edit = str_replace_all(area_edit, '\\s+', ' '))

### step 2: expand 'to' for areas (e.g. areas 3 to 10).  
### * Identify clauses, split into 'from'/'to'
### * create a vector of from:to and convert to string.
### * str_replace the clause with the expanded vector.
### * reassemble the entire closure list string
lic_lookup_to <- lic_lookup_clean %>%
  filter(str_detect(area_edit, '[0-9] to [0-9]')) %>%
  mutate(area_split = str_split(tolower(area_edit), ',|;|and')) %>% ### divide at comma or semicolon
  unnest(area_split) %>%
  mutate(to_clause = str_extract(area_split, '[0-9]*~?[0-9]+ to [0-9]+~?[0-9]*')) %>%
  separate(to_clause, into = c('from_num', 'to_num'), sep = '[^0-9]? to ', remove = FALSE) %>%
  rowwise() %>%
  ### build vector of Area numerals, skipping sub-areas
  mutate(to_vec = ifelse(!is.na(from_num) & !str_detect(from_num, '~') & !str_detect(to_num, '~'),
                         paste(as.integer(from_num):as.integer(to_num), collapse = ', '),
                         '')) %>%
  ### build vector of Sub-Area numerals, skipping areas
  mutate(tmp_area = ifelse(!is.na(from_num) & str_detect(from_num, '~') & str_detect(to_num, '~'),
                           str_extract(from_num, '[0-9]+(?=~)'), NA),
         to_vec_sub = ifelse(!is.na(from_num) & str_detect(from_num, '~') & str_detect(to_num, '~'),
                         paste(tmp_area, as.integer(str_extract(from_num, '(?<=~)[0-9]+')):as.integer(str_extract(to_num, '(?<=~)[0-9]+')), sep = '~', collapse = ', '),
                         '')) %>%
  ungroup() %>%
  mutate(to_vec = ifelse(to_vec == '', to_vec_sub, to_vec),
         area_split2 = ifelse(!is.na(to_clause),
                              str_replace(area_split, to_clause, to_vec),
                              area_split),
         area_split3 = str_replace_all(area_split2, '[^0-9~]+', ' '),
         area_split3 = str_trim(area_split3)) %>%
  select(-area_split, -area_split2, -to_clause, -from_num, -to_num, -to_vec) %>%
  distinct() %>%
  filter(!str_detect(area_split3, '[0-9]{4}')) %>% ### ditches four digit #s (i.e. years)
  group_by(lic_area) %>%
  summarize(desc      = first(desc),
            area_edit = paste(area_split3, collapse = ' ')) %>%
  ungroup()
  
lic_lookup_all <- lic_lookup_clean %>%
  filter(!lic_area %in% lic_lookup_to$lic_area) %>%
  mutate(area_edit = str_replace_all(area_edit, '[0-9]{4}', ''),
         area_edit = str_replace_all(area_edit, '[^0-9~]+', ' '),
         area_edit = str_trim(area_edit)) %>%
  bind_rows(lic_lookup_to) %>%
  mutate(pfma_area = str_split(area_edit, ' ')) %>%
  unnest(pfma_area) %>%
  separate(pfma_area, c('pfma', 'pfmsa'), '~') %>%
  mutate(pfma = as.integer(pfma),
         pfmsa = as.integer(pfmsa)) %>%
  select(-area_edit) %>%
  distinct()

### attach OHIBC regions by PFMA
pfma_lookup <- read_csv(file.path(dir_git, 'prep/fis/v2017/int/pfmsa_to_ohibc.csv')) %>%
  select(rgn_id, pfma_id, pfmsa_id) %>%
  filter(!is.na(pfma_id) & !is.na(rgn_id)) %>%
  distinct()

lic_lookup_rgn <- lic_lookup_all %>%
  filter(is.na(pfmsa)) %>%
  inner_join(pfma_lookup, by = c('pfma' = 'pfma_id')) %>%
  bind_rows(lic_lookup_all %>%
              filter(!is.na(pfmsa)) %>%
              inner_join(pfma_lookup, by = c('pfma' = 'pfma_id', 'pfmsa' = 'pfmsa_id'))) %>%
  select(-pfmsa_id, -pfmsa) %>%
  distinct()

write_csv(lic_lookup_rgn, file.path(dir_goal, 'int/lic_lookup_rgn.csv'))

3.4.1.1 License areas cleaned and expanded

lic_lookup_rgn <- read_csv(file.path(dir_goal, 'int/lic_lookup_rgn.csv')) %>%
  filter(!is.na(rgn_id))

lic_clean <- read_csv(lic_clean_file) %>%
  mutate(area = tolower(area))

lic_fn_by_rgn <- lic_clean %>%
  inner_join(lic_lookup_rgn, by = c('area' = 'lic_area')) %>%
  select(-pfma) %>%
  distinct()

x <- lic_clean %>%
  group_by(area) %>%
  mutate(year_1 = min(year, na.rm = TRUE),
            year_n = max(year, na.rm = TRUE)) %>%
  ungroup()

y <- x %>%
  filter(!year <= 1995) %>%
  filter(!str_detect(area, 'yukon|stikine|taku')) %>%
  filter(!is.na(area))
### still some areas not bound; these all end pre-2003, most pre-1999

sum_fn_by_rgn <- lic_fn_by_rgn %>%
  group_by(year, rgn_id) %>%
  summarize(n_fn = sum(fn),
            n_tot = n(),
            pct_fn = n_fn/n_tot) %>%
  filter(year >= 1995) %>% ### no First Nations licenses before 1995 on record
  filter(year < 2017)      ### incomplete current year

write_csv(sum_fn_by_rgn, file.path(dir_goal, 'output/ao_licenses.csv'))
sum_fn_by_rgn <- read_csv(file.path(dir_goal, 'output/ao_licenses.csv')) %>%
  left_join(get_rgn_names(), by = 'rgn_id') %>%
  mutate(pct_fn = 100 * pct_fn)

plot_fn_by_rgn <- ggplot(sum_fn_by_rgn, aes(x = year, y = pct_fn, group = rgn_name, color = rgn_name)) +
  ggtheme_plot() +
  geom_line(size = 1.5, alpha = 0.7) +
  labs(y = 'Percent of licenses to First Nations',
       color = 'OHIBC Region') +
  scale_color_brewer(palette = 'Dark2')

print(plot_fn_by_rgn)

bc_poly_84 <- sf::read_sf(dsn = dir_spatial %>% path.expand(),
                          layer = 'ohibc_rgn_simple') %>%
  st_transform(4326) %>%
  mutate(rgn_id = as.integer(rgn_id))

sum_fn_by_rgn <- read_csv(file.path(dir_goal, 'output/ao_licenses.csv')) %>%
  mutate(pct_fn = pct_fn * 100) %>%
  filter(year %% 5 == 0)

bc_poly_fn <- bc_poly_84 %>%
  full_join(sum_fn_by_rgn, by = 'rgn_id')

ao_max_pct <- bc_poly_fn$pct_fn %>% max()

bpMap <- ggplot(bc_poly_fn) +
  ggtheme_plot() +
  theme(axis.text = element_blank(),
        axis.title = element_blank()) +
  geom_sf(aes(fill = pct_fn),
          color = 'slateblue', size = .25, alpha = .4) +
  scale_fill_gradientn(colours = brewer.pal(7, 'RdYlBu'), limits = c(0, ao_max_pct)) +
  labs(title="OHIBC AO",
       fill = '% licenses\nFirst Nation') +
  facet_wrap( ~ year)
  
ggsave(filename = file.path(dir_goal, 'int/ao_licenses_mapped.png'),
       plot = bpMap)

3.5 Reference point as % FN in rgn?

One possible reference point is that the % of licenses apportioned to FN should be proportional to the % of population in region that is First Nations.

  • Use census data from 2011 or 2016, mapped using the inland regions (full watersheds).
  • Use population of First Nations communities vs total population.
  • It may be desirable to set a lower limit to ensure some minimal representation of First Nations in a region’s fishing population.
ohibc_rgn <- read_sf(dsn = dir_spatial,
                     layer = 'ohibc_rgns_unclipped')

csd_clean <- function(csd_sf, year) {
  csd_sf2 <- csd_sf %>%
    clean_df_names() %>%
    select(csduid, csdname, csdtype, geometry) %>%
    st_transform(st_crs(ohibc_rgn))
  
  csd_sf2$a_tot <- st_area(csd_sf)
  
  csd_sf_bc <- csd_sf2 %>%
    st_buffer(dist = 0) %>%
    st_intersection(ohibc_rgn) 
  
  csd_sf_bc$a_ohibc <- st_area(csd_sf_bc)
  
  csd_bc_df <- csd_sf_bc %>%
    as.data.frame() %>%
    clean_df_names()  %>%
    mutate(prop_area = as.numeric(a_ohibc / a_tot),
           area_total_km2 = as.numeric(a_tot / 1e6),
           area_ohibc_km2 = as.numeric(a_ohibc / 1e6),
           year = year) %>%
    select(csduid, csdname, csdtype, rgn_id, 
           area_total_km2, area_ohibc_km2, prop_area, year)
  
  return(csd_bc_df)
}

if(!file.exists(file.path(dir_goal, 'int/ao_licenses_csd_2016.csv'))) {
  
  csd_2016 <- read_sf(dsn = file.path(dir_anx, '_raw_data/bcstats/csd2016carto'),
                      layer = 'lcsd000b16a_e') %>%
    csd_clean(2016)
  
  write_csv(csd_2016, file.path(dir_goal, 'int/ao_licenses_csd_2016.csv'))

}
pop_2016 <- read_csv(file.path(dir_anx, '_raw_data/bcstats',
                               '2016 Census - CSD_control_79d4a796-2ccd-4732-a01c-1a8fd26d254a_1.csv')) %>%
  clean_df_names() %>%
  select(sgc, name, type, 
         pop = `2016_population`, 
         prev_pop = `2011_population_2`)

write_csv(pop_2016, file.path(dir_goal, 'int/ao_licenses_pop_2016.csv'))
pop_df <- read_csv(file.path(dir_goal, 'int/ao_licenses_pop_2016.csv'))

csd_df <- read_csv(file.path(dir_goal, 'int/ao_licenses_csd_2016.csv')) 

pop_csd <- pop_df %>%
  left_join(csd_df, by = c('sgc' = 'csduid')) %>%
  filter(!is.na(rgn_id)) %>%
  select(-type, -name, -prev_pop) %>%
  group_by(rgn_id, csdtype) %>%
  summarize(pop = sum(pop * prop_area, na.rm = TRUE))

write_csv(pop_csd, file.path(dir_goal, 'int/ao_licenses_pop_csd_2016.csv'))

3.5.1 Determine population of First Nations communities within OHIBC

Develop a dataframe of municipalities by OHIBC region (inland, based on MaPP inland regions and watershed boundaries). Census subdivisions will be allocated to OHIBC regions according to area-weighting.

First Nations communities can be identified from the census data using the “type” codes:

CSD type Description
CY City
DM District Municipality
T Town
VL Village
IM Island Municipality
IRI ✓ Indian Reserve
RDA Regional District Electoral Area
IGD ✓ Indian Governmental District
S-E ✓ Indian Settlement
NL ✓ Nisga’a Land
NVL ✓ Nisga’a Village

Note, this does not catch First Nations members who are living separately from these communities, and counts all members of the communities as First Nations members.

pop_csd <- read_csv(file.path(dir_goal, 'int/ao_licenses_pop_csd_2016.csv'))

fn_types <- c('IRI', 'IGD', 'S-E', 'NL', 'NVL') %>%
  paste(collapse = '|')

pop_fn_rgn <- pop_csd %>%
  group_by(rgn_id) %>%
  mutate(fn_pop = str_detect(csdtype, fn_types)) %>%
  summarize(fn_pop = sum(pop * fn_pop),
            tot_pop = sum(pop),
            pct_fn_pop = round(fn_pop / tot_pop, 5))


knitr::kable(pop_fn_rgn, caption = 'First Nations communities as percent of total pop')
First Nations communities as percent of total pop
rgn_id fn_pop tot_pop pct_fn_pop
1 3178.346 4.415680e+04 0.07198
2 1435.355 4.328024e+03 0.33164
3 2259.224 3.625213e+03 0.62320
4 2336.893 3.324516e+04 0.07029
5 4937.957 3.328312e+05 0.01484
6 28014.948 3.253829e+06 0.00861
8 0.000 1.009863e+00 0.00000
write_csv(pop_fn_rgn, file.path(dir_goal, 'output/ao_licenses_fn_pop.csv'))
license_ref_min <- 0.20 ### arbitrary lower limit of reference point

license_status <- read_csv(file.path(dir_goal, 'output/ao_licenses.csv')) %>%
  filter(rgn_id != 7) %>%
  left_join(read_csv(file.path(dir_goal, 'output/ao_licenses_fn_pop.csv')), 
            by = 'rgn_id') %>%
  mutate(ref_pt = ifelse(pct_fn_pop > license_ref_min, pct_fn_pop, license_ref_min),
         status = pct_fn / ref_pt,
         status = ifelse(status > 1, 1, status),
         component = 'first_nations_licenses') %>%
  select(year, rgn_id, pct_fn, pct_fn_pop, ref_pt, status, component) %>%
  left_join(get_rgn_names(), by = 'rgn_id')


plot_fn_by_rgn <- ggplot(license_status, aes(x = year, y = pct_fn, group = rgn_name, color = rgn_name)) +
  ggtheme_plot() +
  geom_hline(aes(yintercept = pct_fn_pop), color = 'grey', size = 1.5, alpha = .5) +
  geom_hline(aes(yintercept = ref_pt), color = 'red', alpha = 1) +
  geom_line(size = 1.5, alpha = 0.7) +
  labs(y = 'Percent of licenses to First Nations',
       color = 'OHIBC Region') +
  scale_color_brewer(palette = 'Dark2') +
  facet_wrap( ~ rgn_name)

print(plot_fn_by_rgn)

Note: grey line indicates % of rgn population that are members of First Nations communities according to 2016 census subdivisions; red line indicates reference point as proportion of region population that is First Nations, with a minimum set at 20%.

pop_fn_rgn <- read_csv(file.path(dir_goal, 'int/ao_licenses_fn_pop.csv'))

fn_by_rgn <- read_csv(file.path(dir_goal, 'output/ao_licenses.csv')) %>%
  left_join(get_rgn_names(), by = 'rgn_id') %>%
  left_join(pop_fn_rgn, by = 'rgn_id') %>%
  filter(rgn_id != 7) %>%
  rename(pct_fn_lic = pct_fn)

fn_by_rgn_mean <- fn_by_rgn %>%
  group_by(rgn_id, rgn_name, pct_fn_pop, fn_pop) %>%
  summarize(pct_fn_lic = mean(pct_fn_lic, na.rm = TRUE),
            n_fn = mean(n_fn, na.rm = TRUE)) %>%
  ungroup()

x <- ggplot(fn_by_rgn, aes(x = pct_fn_pop, y = pct_fn_lic)) +
  geom_abline(slope = 1, intercept = 0, color = 'red') +
  geom_point(aes(color = rgn_name), alpha = .3) +
  geom_point(data = fn_by_rgn_mean, aes(color = rgn_name), size = 4) +
  labs(title = 'FN % licenses vs % rgn pop')

y_mdl <- lm(n_fn ~ fn_pop, data = fn_by_rgn)
y_coef <- c('slope' = y_mdl$coefficients[['fn_pop']],
            'intcp' = y_mdl$coefficients[['(Intercept)']],
            'r2'    = y_mdl %>% summary() %>% .$r.squared)

mean_ratio <-  fn_by_rgn %>%
  group_by(rgn_id) %>%
  summarize(mean_n = mean(n_fn),
            mean_fn_pop = mean(fn_pop),
            mean_ratio = mean(n_fn / fn_pop)) %>%
  ungroup() %>%
  filter(rgn_id != 8) %>%
  mutate(unif_dist = sum(mean_n) / sum(mean_fn_pop))

y <- ggplot(fn_by_rgn, aes(x = fn_pop, y = n_fn)) +
  geom_smooth(method = 'lm', color = 'red') +
  geom_point(aes(color = rgn_name), alpha = .3) +
  geom_point(data = fn_by_rgn_mean, aes(color = rgn_name), size = 4) +
  labs(title = 'FN licenses vs rgn pop')

  
print(x); print(y)

3.6 Other questions

3.6.1 licenses by length?

Do licensed vessels differ greatly between the First Nations fleet and the non-First Nations fleet? A significant difference in boat size could indicate a significant difference in the access to catch for one license type vs. another.

boats <- read_csv(lic_clean_file, nogit = TRUE) %>%
  select(year, fishery_group, fn, length = vessel_length_m, lic_holder, lic_oper) %>%
  mutate(length = ifelse(length == 0, NA, length))

bigboats <- boats %>% 
  filter(length > 50)

medboats <- boats %>%
  group_by(year, fn) %>%
  mutate(mean_l = mean(length, na.rm = TRUE)) %>%
  ungroup() %>%
  filter(length <= 40) %>%
  filter(year %in% c(1995, 2000, 2005, 2010, 2015))

length_hist <- ggplot(medboats, aes(x = length, fill = fn)) +
  geom_histogram(aes(frame = year), alpha = .6) +
  geom_vline(aes(xintercept = mean_l, color = fn)) +
  geom_text(data = medboats %>% 
              select(year, mean_l, fn) %>% 
              distinct(),
            aes(x = mean_l, y = 1, label = paste0('mean = ', round(mean_l, 1), ' m')),
            hjust = 0, vjust = 0, angle = 90, size = 3) +
  labs(x = 'vessel length, meters',
       y = 'license count',
       fill = 'First Nations?',
       color = 'First Nations?') +
  facet_wrap( ~ year)

# print(length_hist)

ggsave(file.path(dir_goal, 'int/ao_licenses_vessel_lengths.png'))

3.6.2 licenses by area intensity and type

lic_lookup_rgn <- read_csv(file.path(dir_goal, 'int/lic_lookup_rgn.csv')) %>%
  filter(!is.na(rgn_id))

lic_clean <- read_csv(file.path(dir_goal_anx, 'int/licenses_clean.csv')) %>%
  mutate(area = tolower(area))

lic_by_pfma_all <- lic_clean %>%
  inner_join(lic_lookup_rgn, by = c('area' = 'lic_area')) %>%
  select(-rgn_id) %>%
  distinct()

### Total licenses by fishery group:
# lic_by_pfma_all$fishery_group %>% table()
#             CLAM                    CRAB  GEODUCK AND HORSE CLAM                 HERRING   HERRING - FOOD & BAIT HERRING - SPECIAL ISSUE 
#           113384                   48600                    4024                      74                    4760                     348 
# OYSTERS, PACIFIC              PROCESSING          RED SEA URCHIN                ROCKFISH    ROE HERRING GILL NET       ROE HERRING SEINE 
#             1650                     283                   14776                  331645                   89560                   16515 
#  SALMON GILL NET       SALMON HISTORICAL            SALMON SEINE            SALMON TROLL            SEA CUCUMBER 
#           316457                      20                   97852                  208194                    6047

### FN licenses by fishery group:
# lic_by_pfma_all %>% filter(fn) %>% .$fishery_group %>% table()
#              CLAM                   CRAB GEODUCK AND HORSE CLAM       OYSTERS, PACIFIC         RED SEA URCHIN               ROCKFISH   ROE HERRING GILL NET 
#             24405                   1898                      9                    156                   1266                  12061                   3750 
# ROE HERRING SEINE        SALMON GILL NET           SALMON SEINE           SALMON TROLL           SEA CUCUMBER 
#               236                  75222                   4271                   6914                     10 

### These have NO fn licenses:
# lic_by_pfma_all %>% group_by(fishery_group) %>% filter(sum(fn) == 0) %>% .$fishery_group %>% table()
# HERRING   HERRING - FOOD & BAIT HERRING - SPECIAL ISSUE              PROCESSING       SALMON HISTORICAL 
#      74                    4760                     348                     283                      20 

lic_by_pfma <- lic_by_pfma_all %>%
  select(year, lic_tab, fn, fishery_group, desc, pfma) %>%
  mutate(taxa = case_when(str_detect(fishery_group, 'CLAM|OYSTER') ~ 'bivalve',
                          str_detect(fishery_group, 'SALMON')      ~ 'salmon',
                          str_detect(fishery_group, 'HERRING')     ~ 'herring',
                          str_detect(fishery_group, 'CRAB')        ~ 'crab',
                          str_detect(fishery_group, 'URCH|CUCUM')  ~ 'echinoderm',
                          str_detect(fishery_group, 'ROCKFISH')    ~ 'rockfish',
                          str_detect(fishery_group, 'PROCESS')     ~ 'processing',
                          TRUE                                     ~ 'other')) %>%
  group_by(year, pfma, taxa, fn) %>%
  summarize(n_lic = n()) %>%
  ungroup()

write_csv(lic_by_pfma, file.path(dir_goal, 'int/lic_by_pfma.csv'))

lic_by_year <- lic_by_pfma %>%
  group_by(year, taxa) %>%
  summarize(fn_lic = sum(n_lic * fn),
            tot_lic = sum(n_lic)) %>%
  ungroup() %>%
  filter(year >= 2000 & year <= 2016) %>%
  filter(taxa != 'processing')

plot_by_year <- ggplot(lic_by_year, aes(x = year)) +
  ggtheme_plot() +
  geom_line(aes(y = fn_lic, group = taxa), 
            color = 'darkred', alpha = .7, size = 1) +
  geom_line(aes(y = tot_lic, group = taxa), 
            color = 'darkblue', alpha = .5, size = 1) +
  scale_color_brewer(palette = 'Dark2') +
  facet_wrap( ~ taxa, scales = 'free_y') +
  ylim(c(0, NA))

print(plot_by_year)

ggsave(file.path(dir_goal, 'figs/licenses/lic_by_year_and_taxa.png'))
lic_by_pfma_tot <- read_csv(file.path(dir_goal, 'int/lic_by_pfma.csv')) %>%
  filter(year >= 2000 & year <= 2016) %>%
  group_by(pfma, taxa) %>%
  summarize(fn_lic = sum(n_lic * fn),
            tot_lic = sum(n_lic)) %>%
  ungroup()


pfma_sf <- sf::read_sf(dsn = file.path(dir_anx, '_raw_data/dfo_khunter', 
                                 'management_boundaries/d2016',
                                 'pac_fishery_mgmt_areas'),
                       layer = 'DFO_BC_PFMA_AREAS_50K_V3_1') %>%
  select(pfma = STATAREA)

bbx <- st_bbox(pfma_sf)[c(1, 3, 2, 4)] %>%
  raster::extent()

garbage <- capture.output({
  continent_sf <- readOGR(dsn = path.expand(file.path(dir_spatial)),
                          layer = 'ohibc_continent') %>%
    raster::crop(bbx) %>%
    sf::st_as_sf()
})


pfma_lic_sf <- pfma_sf %>%
  left_join(lic_by_pfma_tot, by = 'pfma') %>%
  filter(!is.na(taxa))

for(taxa_sel in unique(pfma_lic_sf$taxa)) {
  ### taxa_sel <- 'herring'
  # cat(taxa_sel, '\n')
  
  tmp_df <- pfma_lic_sf %>%
    filter(taxa == taxa_sel) %>%
    mutate(tot_lic = ifelse(tot_lic == 0, NA, tot_lic),
           fn_lic = ifelse(fn_lic == 0, NA, fn_lic))
  
  breaks <- c(1, 10, 100, 1000, 10000)
  max_tot <- tmp_df$tot_lic %>% max(na.rm = TRUE)
  breaks_tot <- breaks[breaks < max_tot * 10]
  
  tot_map <- ggplot(tmp_df %>% filter(!is.na(tot_lic))) +
    ggtheme_plot() +
    geom_sf(data = continent_sf, fill = '#bb9999', color = NA) +
    geom_sf(data = pfma_sf, fill = 'grey96', 
            color = 'black', size = .1) +
    geom_sf(aes(fill = tot_lic), 
            color = NA, alpha = .8) +
    scale_fill_distiller(palette = 'Blues', 
                         breaks = breaks_tot,
                         limits = c(1, NA),
                         trans = 'log', direction = 1) +
    labs(title = taxa_sel,
         fill  = 'Total licenses')
  
  print(tot_map)

  ggsave(file.path(dir_goal, sprintf('figs/licenses/lic_map_pfma_%s_total.png', taxa_sel)))

  tmp_df_fn <- tmp_df %>% filter(!is.na(fn_lic))
  if(nrow(tmp_df_fn) > 0) {
    max_fn  <- tmp_df$fn_lic  %>% max(na.rm = TRUE)
    breaks_fn  <- breaks[breaks < max_fn  * 10]

    fn_map <- ggplot(tmp_df_fn) +
      ggtheme_plot() +
      geom_sf(data = continent_sf, fill = '#bb9999', color = NA) +
      geom_sf(data = pfma_sf, fill = 'grey96', 
              color = 'black', size = .1) +
      geom_sf(aes(fill = fn_lic), 
              color = NA, alpha = .8) +
      scale_fill_distiller(palette = 'Blues', 
                           breaks = breaks_fn,
                           limits = c(1, NA),
                           trans = 'log', direction = 1) +
      labs(title = taxa_sel,
           fill  = 'FN licenses')
  
    print(fn_map)
    ggsave(file.path(dir_goal, sprintf('figs/licenses/lic_map_pfma_%s_fn.png', taxa_sel)))
  }
}


prov_wrapup(commit_outputs = FALSE)

4 Provenance

  • Run ID: 15 (30ed76d); run tag: “standard run”
  • Run elapsed time: 988.369 seconds; run memory usage: 61725.5 MB
  • System info:
    • System: Linux, Release: 4.4.0-57-generic. Machine: x86_64. User: ohara.
    • R version: R version 3.4.1 (2017-06-30), Platform: x86_64-pc-linux-gnu, Running under: Ubuntu 14.04.5 LTS.
    • Attached base packages: stats, graphics, grDevices, utils, datasets, methods, base
    • Other attached packages: bindrcpp_0.2, provRmd_0.1.1, stringr_1.2.0, RColorBrewer_1.1-2, dplyr_0.7.2, purrr_0.2.2, readr_1.1.0, tidyr_0.6.3, tibble_1.3.3, ggplot2_2.2.1.9000, tidyverse_1.1.1, broom_0.4.2, sf_0.5-1, rgdal_1.2-8, sp_1.2-5
%3 16 Pacific Region Commercial Fishing Licenses En.csv 12 ao_licenses_prep.rmd#read_xlsx 16->12 used 17 Commercial License Meta Data EN (1).csv 17->12 used 18 licenses_clean.csv 4 ao_licenses_prep.rmd#collate_fishery_groups 18->4 used 13 ao_licenses_prep.rmd#total_licenses 18->13 used 14 ao_licenses_prep.rmd#unnamed-chunk-1 18->14 used 3 ao_licenses_prep.rmd#check_areas_vs_licenses 18->3 used 5 ao_licenses_prep.rmd#examine_license_area_intensity 18->5 used 19 license_area_lookup.csv 19->14 used 20 pfmsa_to_ohibc.csv 20->14 used 21 lic_lookup_rgn.csv 21->3 used 21->5 used 22 ao_licenses.csv 8 ao_licenses_prep.rmd#plot_fn_by_rgn 22->8 used 15 ao_licenses_prep.rmd#use_ggplot_and_sf 22->15 used 9 ao_licenses_prep.rmd#plot_vs_pop_reference 22->9 used 23 ao_licenses_mapped.png 24 2016 Census - CSD_control_79d4a796-2ccd-4732-a01c-1a8fd26d254a_1.csv 11 ao_licenses_prep.rmd#read_pop_data 24->11 used 25 ao_licenses_pop_2016.csv 2 ao_licenses_prep.rmd#census_pops_2016 25->2 used 26 ao_licenses_csd_2016.csv 26->2 used 27 ao_licenses_pop_csd_2016.csv 7 ao_licenses_prep.rmd#ohibc_pop_fns 27->7 used 28 ao_licenses_fn_pop.csv 28->9 used 29 lic_by_pfma.csv 6 ao_licenses_prep.rmd#map_lic_by_pfma_and_year 29->6 used 30 lic_by_year_and_taxa.png 31 ohibc_continent.shp 31->6 used 32 lic_map_pfma_bivalve_total.png 33 lic_map_pfma_bivalve_fn.png 34 lic_map_pfma_crab_total.png 35 lic_map_pfma_crab_fn.png 36 lic_map_pfma_echinoderm_total.png 37 lic_map_pfma_echinoderm_fn.png 38 lic_map_pfma_rockfish_total.png 39 lic_map_pfma_rockfish_fn.png 40 lic_map_pfma_salmon_total.png 41 lic_map_pfma_salmon_fn.png 42 lic_map_pfma_herring_total.png 43 lic_map_pfma_herring_fn.png 1 ao_licenses_prep.rmd 1->12 wasExecutedBy 1->4 wasExecutedBy 1->13 wasExecutedBy 1->14 wasExecutedBy 1->3 wasExecutedBy 1->8 wasExecutedBy 1->15 wasExecutedBy 1->11 wasExecutedBy 1->2 wasExecutedBy 1->7 wasExecutedBy 1->9 wasExecutedBy 1->5 wasExecutedBy 1->6 wasExecutedBy 10 ao_licenses_prep.rmd#prov_footer 1->10 wasExecutedBy 12->18 wasGeneratedBy 14->21 wasGeneratedBy 3->22 wasGeneratedBy 15->23 wasGeneratedBy 11->25 wasGeneratedBy 2->27 wasGeneratedBy 7->28 wasGeneratedBy 5->29 wasGeneratedBy 5->30 wasGeneratedBy 6->32 wasGeneratedBy 6->33 wasGeneratedBy 6->34 wasGeneratedBy 6->35 wasGeneratedBy 6->36 wasGeneratedBy 6->37 wasGeneratedBy 6->38 wasGeneratedBy 6->39 wasGeneratedBy 6->40 wasGeneratedBy 6->41 wasGeneratedBy 6->42 wasGeneratedBy 6->43 wasGeneratedBy